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Abstract 

In this article we prove that it is possible to construct, using newest- 
vertex bisection, meshes that equidistribute the error in _ff 1 -norm, when- 
ever the function to approximate can be decomposed as a sum of a reg- 
ular part plus a singular part with singularities around a finite number 
of points. This decomposition is usual in regularity results of Partial Dif- 
ferential Equations (PDE). As a consequence, the meshes turn out to be 
quasi-optimal, and convergence rates for adaptive finite element meth- 
ods (AFEM) using Lagrange finite elements of any polynomial degree are 
obtained. 

1 Introduction 

Adaptive procedures for the numerical solution of partial differential equations 
(PDE) started in the late 70's and are now standard tools in science and en- 
gineering. The ultimate purpose of adaptivity is to reduce the computational 
cost through the automatic construction of a sequence of meshes that would 
eventually equidistribute the approximation errors, leading to (quasi-)optimal 
meshes. Adaptive methods for stationary problems usually consist of the loop 

SOLVE -> ESTIMATE -> MARK -► REFINE. (1) 

Experience strongly suggests that, starting from a coarse mesh, such an 
iteration converges within any prescribed error tolerance in a finite number of 
steps, and it does so in an optimal manner, provided the a posteriori error 
estimators are reliable and efficient. What is observed in fact, is that for a 
large class of problems and data, the solutions uq- and meshes T obtained with 
adaptive methods of the form (1) satisfy 

\\u-ur\\m <C(#T)-^, (2) 



where u denotes the exact solution, p the polynomial degree of the finite el- 
ement space over the mesh T, and d the dimension of the underlying space. 
This is the same error bound that is obtained with uniformly refined meshes 
for smooth (regular) solutions u G H p+1 , by an application of classical in- 
terpolation estimates [Ciarlet 1978]. The decay rate dictated by (2) — which 
is also observed in practice for the so-called singular solutions belonging to 
H s (£l) for s < 2 — is usually called optimal error decay. The precise goal of 
this paper is to show a broader family of functions for which this so-called 
optimal decay can be obtained when using adaptive methods. We will prove 
that this decay holds for functions that can be decomposed as a sum of a reg- 
ular part plus singular terms, as described in classical regularity results for 
PDE [Grisvard 1985, Grisvard 1992, Petzoldt 2001, Dauge 1988]. 

The first steps towards understanding the optimality of AFEM consisted 
of studying their convergence. An analysis of (1) for linear, elliptic, and sym- 
metric problems in Id is presented in [Babuska Vogelius 1984]. The first mul- 
tidimensional result is given in [Dorflcr 1996], where it is proved that, after 
a pre-adaptation to data, (1) reduces the error below any prescribed toler- 
ance. Proper convergence without conditions on the initial grid is proved 
in [Morin Nochctto Siebcrt 2002], requiring the so-called interior node prop- 
erty and an additional marking step driven by data oscillation. The latter 
work was generalized in various directions. Lately, convergence of adaptive 
methods with marking strategies other than Dorfler's, for a large class of lin- 
ear problems with different a posteriori error estimators, and without requir- 
ing the marking due to oscillation or the interior node property, was proved 
in [Morin Siebert Veeser 2007]. The result only leads to asymptotic conver- 
gence without an error reduction in every step, which seems to be essential to 
prove optimality though (see [Stevenson 2006, Cascon ct. al. 2007]). 

Regarding complexity, an important result for an algorithm which is very 
similar to (1), is proved in [Stevenson 2006]. The proof relies on techniques first 
developed in [Bincv Dahmen DeVore 2004] and new ideas. This result was later 
improved in several aspects in [Cascon et. al. 2007]: the artificial assumptions of 
interior node and marking due to data oscillation were removed, and the result 
applies to more general elliptic equations. 

When considering adaptive methods the notion of complexity differs from 
the previous one which was based on a uniform clement size h. It is now defined 
in terms of the number of elements (or degrees of freedom) necessary to achieve 
a certain tolerance. 

In order to be more specific at this point we need to introduce some notation. 
Let us assume that we have a function u £ H 1 ^), where 17 is a polygonal 
domain in M 2 (polyhedral in R 3 ), and H 1 (f2) denotes the Sobolev space of square 
integrable functions with square integrable weak derivatives of first order. 

We consider an initial triangulation Tq of the domain Q into simplices, and we 
let the admissible triangulations be those obtained from Tq with newest- vertex 
bisection, either the iterative [Bansch 1991] or the recursive [Kossaczky 1994] 
version, without hanging nodes. For each admissible triangulation T we consider 
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the Lagrange finite element space 



V r = {v G H\n) : v\ T G P P ,VT e T} , 

where, for p G N, V p denotes the space of polynomials of degree < p. The best 
approximation error with complexity N, for N G N, is defined as follows: 

a p N {u)= min inf \\u-v\\ H i m , 

where TV := {T admissible : (#T — #7o) < iV} that is, the minimum over T 
is taken over all admissible triangulations obtained with at most N bisections. 
We now define, for s > the approximation classes 

A* = {veH 1 (Sl): 3C such that a p N {v) < CN~ S ,VN G N} , 

or, equivalcntly, 

A p s = {v G : |w| A p < 00} with |v| a p := sup cr p N (v)N s . 

The first complexity results for adaptive finite clement methods (AFEM) are 
presented in [Binev Dahmcn DeVore 2004], for an algorithm that needs coars- 
ening, which seems not to be necessary, at least for symmetric elliptic problems. 
This, and the aforementioned papers on optimality of AFEM [Stevenson 2006, 
Cascon et. al. 2007] study adaptive algorithms for approximating the solution u 
to an elliptic partial differential equation. Essentially, the following fundamen- 
tal result is proved: the adaptive algorithms generate a sequence {(Tk,Uk)}k&i 
of triangulations and finite element approximations uu G V^- fc that satisfy the 
following: 

If u G Al then ||u - u k \\ m{n) < C{#T k )~ s , Vfc G N. 

That is, the sequence of triangulations and approximate solutions have a com- 
plexity with the same decay rate as the optimal ones. The interesting aspect 
of those results is the fact that such a (quasi-)optimal approximation is ob- 
tained through a standard adaptive loop for the elliptic problem, without a 
priori knowledge of the exact solution, and with a number of operations propor- 
tional to the cardinality of the meshes. Notice that a simple minded approach 
to compute <J p N {u) with precise knowledge of u could lead to exponential work 
in terms of N . 

The question — already raised in [Cascon et. al. 2007] — that is still unan- 
swered is what rate s is to be expected in different situations. From the results 
just described it is clear that AFEM do a quasi-optimal job among all possible 
adaptive meshes. What we present in this article, is quantitative information 
about the convergence rate of AFEM. In order to do so, we relate the mem- 
bership of a function to an approximation class A p with its regularity, proving 
rigorously, through the construction of specific meshes, that certain class of 
functions is contained in A p . 
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In [Binev et. al. 2002] an almost characterization of these classes is obtained, 
for the case p = 1 in terms of Besov regularity for Lipschitz polygonal do- 
mains; the proof is based on an adaptive tree approximation algorithm. To 
illustrate the applicability of this result we just mention without giving too 
much detail — that the Besov space -B 2 (L T (f2)) is contained in A]y 2 for all 
r > 1 [Binev ct. al. 2002, Theorem 5.1]. 

The regularity of solutions to Poisson's problem on Lipschitz domains, in 
terms of Besov regularity is studied in [Dahlke DeVore 1997]. It is proved that 
for Poisson's equation —Au = f in a Lipschitz polygonal domain del 2 , with 
homogeneous Dirichlet boundary values, u G _B 2 (L r (£l)) if / £ 

Combining these two results we obtain that u € A]y 2 if / G H 1 , but a 
stronger result holds. Using Grisvard's Sobolev regularity results [Grisvard 1985] , 
we have that only assuming / G L 2 (f2), u G W3(f2), that is, all derivatives of 
order up to two are in L p (tt), for all 1 < p < 4/3. This, in turn implies that 
for all 1 < r < 4/3, u belongs to the Besov space i? 2 (L T (f2)), and applying the 
result [Binev et. al. 2002] this implies u G A]y 2 under the sole assumption of 

/ G L 2 (ft). 

The spirit of the results that we present in this article is a combination 
of [Binev et. al. 2002] and [Dahlke DeVore 1997]. However, our approach will 
not hinge upon regularity in Besov terms but rather upon a decomposition of the 
functions as a sum of a regular part plus singular terms, as stems from the classi- 
cal regularity results for PDE like those stated in [Grisvard 1985, Grisvard 1992, 
Kellogg 1975, Kellogg 1992, Petzoldt 2001, Dauge 1988]. We obtain results for 
polygonal domains which are not necessarily Lipschitz (including slit domains) 
and we generalize to any polynomial degree p > 1; the proof is elementary, and 
does not make use of sophisticated theory of L q spaces for q < 1, as seems nec- 
essary in the approach of [Binev et. al. 2002]. Moreover, our result is directly 
applicable in some cases where the Besov regularity of the solutions to the PDE 
is not available, but instead, a descomposition into a regular plus a singular 
part is known to hold. 

In [Grisvard 1985, Grisvard 1992] one can find some conditions on the cl- 
ement sizes relative to the distance to the points where the singularities are 
located, in order to obtain an error of order iV -1 / 2 when using linear elements 
in 2d. The difference between our result and those, is that we present an al- 
gorithm for constructing those meses using bisection, and thus show that those 
meshes are attainable by an adaptive algorithm. Moreover, in view of the results 
in [Stevenson 2006, Cascon et. al. 2007], a consequence of our result is that the 
standard adaptive algorithms proposed there generate a sequence of meshes and 
discrete solutions {7k,Uk} k satisfying \\u — UkWn 1 ^) < C(#^fc) • A quanti- 
tative answer regarding convergence rates of adaptive finite element methods is 
thus obtained, for Lagrange finite elements of any polynomial degree p > 1. 

The rest of the article is organized as follows. In section 2 we state the main 
result and present some applications to solutions of elliptic PDE in section 3. In 
section 4 we propose an algorithm for constructing the desired mesh and prove 
some of its properties. We conclude the proof of the main result by bounding 
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the error in section 5. 



2 Main Result 

From now on, for any admissible triangulation T of the domain £1, we let Vt 
denote the finite element space of continuous piecewise polynomials of degree 
< p, where p is a fixed positive integer. The following is the main result of this 
article, which states that a large family of functions, as those obtained when 
solving elliptic and other PDE, belong to &? p / 2 - 

Theorem 2.1. Let Q C M. d be a polygonal (d = 2) or polyhedral (d = 3) domain, 
not necessarily Lipschitz, let Tq be an initial triangulation of£l and suppose that 

N 

u = ^ (3) 

i=Q 

where: 

• u a G H^fl), with u \t G HP +1 {T), for all T G T Q ; 

• for i = 1, 2, . . . , TV, Ui can be expressed in polar coordinates around Xi as 

ui = Q (ln(r;)) fc * r]' gt(di)xi, 

where: 

1. Cj are constants and ki are nonnegative integers. 

{ x i}iLi —'■ N is a set of points in f2, that are also vertices of To; 
3. Ti denotes the distance to a;,-, and: 

— 9i = 8i G [0, 2n) is the angle coordinate of x with respect to Xi 
and a half line starting at Xj, when d = 2; 

— 6i = (6i,(f>i) G [0, 2n) x [0, 7r], where <j>i is the angle coordinate of 
x with respect to Xi and a half line R starting at Xi , and letting 
P denote the plane orthogonal to R that contains the 
angle coordinate of the projection of x on the plane a half line S 
starting at X{ contained into P, when d = 3. 

4- 7i are positive constants; 

5. the functions g% satisfy the following assumptions depending on the 
dimension d: 

— g% £ W^(0, 2ir), satisfies the periodicity condition <?i(0) = gi(2ir) 
and is piecewise W^ 1 in the following sense: there exists a par- 
tition of [0, 2ir] into segments such that gi\s G W^~ 1 (S) for 
all S G tyi, when d = 2; 
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— gi G W^q^O, 2ir) x (0, n)), satisfies the periodicity conditions 
ffi(O,0i) = 5i(27r, < & < 7r, and ffi(0,0) = gi(0i,0), 
<7i(0, 27r) = gi(0i,2ir), < 0, < 27r, and is piecewise W^ 1 in 
the following sense: there exists a partition tyi o/ (0, 27r) x (0, 7r) 
into triangles such that gi\s G W^' 1 (S) for all S G tyi, w/ien 
d = 3; 

6. x« are C°°(Q) cutoff functions; 

7. the jumps of Vu^ (if any) are aligned with the edges of the initial 
mesh Tq. 

Then, for any given tolerance e > 0, there exists a conforming triangulation T , 
obtained by newest-vertex bisection, starting from 7q such that: 
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inf ||u-«r|| in <e and # r - # r o < C U , T „ — , (4) 



where C u ,t depends on all the parameters that enter the definition of the singu- 

and 

1/2 



lar part X)i=i u i> 071 ^0> arl ^ on u through the broken seminorm \uq\ h p^ 



\ 1 / z 

StgTo ll-^' P+1 ' lt o|||2(x) J , but not on e. Therefore u G ^p/d- 

It is worth observing that, if u satisfies the assumptions of the theorem, then 
we can only assure that u G H 1+e (fl) for all < e < minx<i<jv 7i- Uniform 
global refinements would only lead to u G A^y d , but e could be very small, and 
this rate is very pessimistic with respect to the one that can be obtained with 
adaptivity. 

Remark 2.2. In order to shed some light on the assumptions of the theorem, 
we notice that they imply the following: 

• If we let 7 = min j li , we are able to control the singular terms through the 
following bound, 

Cr7>ln(r 4 )V- (5) 

and similar ones. They imply that, for each of the singular terms itj, 
i = 1, 2, . . . , N, there exists a constant C , such that 

N<CY7, |Vm| <Cr]~\ and \DP +1 Ui \ < Cr7~ p ~\ (6) 

the last inequality holding only in the interior of the elements of Tq , and 
thus also in the interior of any element of any refinement of Tq. The 
constant C depends on a, ki, ji, the W / ^+ 1 -norm of Xi, the W^-norm of 
gi, and the piecewise W^ +1 -norm of gi, that is, on the W^ 1 (S')-norm of 
gi, for all S G %. 

The factor | in the definition of 7 is imposed to control the logarithmic 
terms. If all hi = 0, i = 1, . . . , N, then 7 could be chosen equal to mim 7,, 
and the same bounds would hold. 
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• if T is any refinement of To, and T G T with T n A/" = then Uj|t G 
fff+ 1 (T),* = 0,l,...,7V; 



• since p > 1, and c? < 3, the Sobolcv embedding theorem and the fact 
that 7j > i = 1, 2, . . . , N imply that each component U{, i = 0, N, is 
continuous in 17, and consequently also u is continuous; 

This consequences of the assumptions are the main ingredients that will be used 
in the proof of our results below. 

Notation 2.3. From now on, the letter C will denote a constant, not al- 
ways equal, depending on the given function u of the assumption of theo- 
rem 2.1, through the H 1 ^) -norm of uq, the broken seminoma \uo\ H P+i, n ^ := 



gular terms m, i = 1, 2, . . . , N of u as in the second item of the previous remark. 
We will reserve the notation a < b to denote a < cb with a constant c depending 
only on shape regularity, or the geometry of the domain. And a ~ b will indicate 
that a < b and b < a. 

3 Applications 

In this section we state two applications to elliptic PDE in two dimensions in 
order to illustrate the applicability of our result. 

3.1 Poisson Equation 

Let O be a polygonal domain in K 2 , not necessarily Lipschitz. And let u be the 
(weak) solution to 



As a consequence of Theorem 3.1 in [Kellogg 1992] (see also [Dauge 1988], or 
Thm. 3.1 in [Nochctto Vccscr Vcrani 2007]) it holds that if / G H p - 1+e (n) for 
some e > 0, then u can be written as in the assumptions in theorem 2.1, where 
Af = {xi}f =1 is the set of vertices of fi, and fc, = 0, i = 1, 2, . . . , N. 

In the case of p = 1, e can be taken to be zero, i.e. / G L 2 (f2), the set Af 
contains only the vertices of O with inner angle u>i greater than tt (c; = for 
the other vertices), and gi(t) = sm(Trt/u>i) for all i = 1, 2, . . . , N. 

In the case of p > 1, the set Af contains all the vertices of il. In order 
to avoid the pathological cases where at least one inner angle a of O satisfies 
ap/ir G N, we assume that / G H p ~ 1+e (£l) for some e > instead of i/ p_1 (il), 
but this is not such a big restriction in practice. Moreover, this hypothesis can 
be weakened and ask that / G L 2 {£1) and f\ T G H p ~ 1+e {T), for all T G %. 

We conclude then that if / G i/ p_1+e (f2) (piecewise over Tq) then the solu- 
tion u to Poisson's equation (7) belongs to A p , 2 . 
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and the parameters and functions defining the sin- 



Am = /, in SI, 
u = 0, on dil 



(7) 
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3.2 Interface Problems for the Laplacian 



Let fl be a polygonal domain, not necessarily Lipschitz, that can be decomposed 
into disjoint subdomains f2,-, i = l,...,Ud with polygonal boundaries: Q = 
U^Jk. We define the interface T = (U^dJl* \ <9ST)). 

Denote with a(x) = Yl7=i a «XO; {x) the global weight function, which is 
constant and positive on each subdomain fij. 

We want to solve the following problem written in variational form: 

Find u € V : / aVu • Vv dx = / fv dx, Wv E V, (8) 
Jn Jn 

where / € L 2 (fl), V = H^Q) = {v E H^fl) : v\r D = 0}, T D C dfl is the 
Dirichlet boundary. This problem is usually called the interface problem for the 
Laplacian and corresponds to the following strong form 

/, in fi;, i = 1,2, . . . ,n d , 
0, on T D 

0, onTiv =dn\T D , 
8u\q . 

— a,j— - — 1 on dCli PI dCli , 
dn 3 

where n denotes the outer unit normal to f2, and rij that of Qi. 

Following the original ideas from [Kellogg 1992], Pctzoldt proved (see Chap- 
ter 2 in [Petzoldt 2001] and references therein) that the solution u to (8) satisfies 
the assumptions of theorem 2.1 for p = 1, if the mesh Tq matches the bound- 
aries of the subdomains Qi and the points on d£l where the boundary condition 
changes are vertices of To . The points xi correspond to the vertices of the inter- 
face r, <9f2, and to those points on <90 where the boundary condition changes. 

We conclude that if / 6 L 2 {Vl) then by theorem 2.1 the solution u belongs 
to Ajy 2 , and the optimal error decay is recovered. 

It is worth mentioning that for certain singular points xg, the value of 
can be as close to zero as desired, depending on the values of a(x) around xg, 
providing very singular examples for the classical theory. In order to illustrate on 
this, we writedown the formulas derived by Kellogg [Kellogg 1975] to construct 
an exact solution of an elliptic problem with piecewise constant coefficients and 
vanishing right-hand side /; for the particular case SI = ( — 1, l) 2 , a = a\ in the 
first and third quadrants, and a = 02 in the second and fourth quadrants. An 
exact solution u to (8) for / = (and non-homogeneous Dirichlet boundary- 
values) is given in polar coordinates by u(r, 0) = r 7 /z(0), where 

if < 9 < tt/2 

if 7T/2 < 6 < 7T 

if 7T < 9 < 3tt/2 
if 3tt/2 < 6 < 2tt 



-V • (a(x)Vu) 



11 
du 
dn 

dm 



cos((tt/2 - (7)7) ■ cos((0 - tt/2 + p)j) 
cos(pj) ■ cos((# — 7r + (7)7) 

C0S((77) • C0S((# — 7T — p)j) 

cos((tt/2 - 0)7) • cos((e - 37i72 - (7)7) 
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and the numbers 7, p, a satisfy the nonlinear relations 



' R := ai/a 2 = — tan((vr/2 — (7)7) ■ cot(/ry) 
1/R = — tan(p7) • cot(cr7) 
R = — tan(<77) • cot((7r/2 — p)j) 
< 7 < 2 

max{0, 7T7 — 7r} < 2jp < min{7T7, n} 
max{0, 7r — 7T7} < —2717 < min{7r, 2n — 717}. 



(9) 



Choosing 7 = 0.1, and solving (9) for R, p and a using Newton's method we 
obtain R = a 1 /a 2 = 161.4476, p = 7r/4, a = -14.92256. A smaller 7 would lead 
to a larger ratio i?, but in principle 7 may be as close to as desired. 

This function u belongs to the Sobolev space H 1+1 (fl), and is thus barely 
in but — according to our results — still in for all p > 1. That is, 

an adaptive finite clement approximation to a solution like this, using Lagrange 
finite elements of degree p will lead to a sequence of meshes and discrete solutions 
{T k ,u k } k satisfying \\u - Uk\\H l (n) < C (#T k )~ p/2 . On the other hand, the 
Besov regularity of the solutions to (8) is not well established, and thus the 
results of [Binev et. al. 2002] are not yet applicable to the interface problem 
for the laplacian. Until the Besov regularity of solutions to PDE is further 
developed, our result — which is far from being a near characterization of the 
class of functions that can be approximated with optimal decay N~d — still 
provides a useful tool to investigate the convergence rate of AFEM for PDE. 



4 Construction 

From now on we assume that u is as in the assumptions of theorem 2.1 and 
we will present an algorithm to construct via newest-vertex bisection a mesh 
fulfilling the properties stated in the theorem. 

Before we introduce the algorithm we will present a heuristic idea with the 
ideal properties that the optimal mesh should have. This will motivate the 
precise definition of the algorithm, which is rather technical, but achieves with 
controlled complexity the goal of equidistribution. 

4.1 Heuristic Idea 

Everything in this subsection will be heuristics, and is presented here — following 
the arguments in [Grisvard 1985, Liao Nochetto 2002, Babuska et. al. 1996]— 
in order to motivate the properties that the optimal mesh should fulfill. The 
precise, rigorous proof will be given in the following sections, after the algorithm 
for constructing the mesh has been presented. 

In order to introduce the basic idea consider the simplest case of a function u 
written in polar coordinates as u = r 1 sin(7#) on a two dimensional domain with 
a reentrant corner of inner angle tt/j at the origin. Suppose that we approximate 
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u with continuous piecewise linear finite elements (p = 1) on a mesh T. The 
H 1 seminorm \u — Itu\\_t of the error between u and its Lagrange interpolant 
Itu on each element is bounded by h\\D 2 u\\ L 2(T) if T and by ||Z?u|| £ 2( T ) if 
6 T. These quantities (squared) also satisfy the following: 

h 2 T \\D 2 u\\ 2 LHT) 44 (7 - 2) |T| S 44 (7 - 2) , if ^ T, 

H^MIl^T) = / T r 2{l ~ l) rdr S ft 27 , if e T, 

Jo 

where r-r denotes the distance of T to the origin and h T := Tj 1 / 2 ^ diam(T). 
In order to achieve the equidistribution of the local error bounds we then require 
for the mesh T that, given a small parameter h > 0, the elements satisfy 

?44 (7 " 2) ft 27 , if £ T, and ft T ^ ft, if £ T. 

Suppose now that this goal is achievable. More precisely, we can classify the 
elements into rings at dyadic distance to the origin, by defining 

D k = {T eT : 2- k - 1 <r T < 2~ fc } , 

for k 6 N, k < K := |log 2 (l/ft)J, and D K = {T e T : r T < 2~ K ). 

Then, on the one hand, the elements T <E Dk, have size \T\ = = 
/!7 r -(T- 2 ) W2 k ^- 2 \ and thus #L> fe = ^ 2*^-2) = /i" 7 2" fe7 which implies 
that 

#r = 2 #^ fc = h-~< 2- fc7 s h~\ 

k<K k 

On the other hand, the error satisfies 

\u - u h \l n - <* ft^ft 27 = ft 7 - (#T)- 1 . 

And this finally implies that \u — Uh\i,n % (WT) ^ \ an d thus u £ A^y 2 . 

In the case <i = 3 if u has a singularity like r 7 as in the previous example, 
the bound \u — u^li.n < (t^^~) > would be obtained if 

/44 (7_2) = ft 27+1 , if £ T, and /i r ^ ft, if £ T. 
4.2 Algorithm 

In this section we will introduce the algorithm that will achieve using newest- 
vertex bisection, a mesh with the precise grading stated in the previous subsec- 
tion, generalized to polynomials of degree p. 
From now on we will use the notation 

rx = min dist(a;j,X) 
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defined for X compact, typically X is a triangle T or a point x, where N 
denotes the finite set where the singularities are located (as in the assumptions 
of theorem 2.1). 

We choose and fix 7 = min j 7i . This choice allows us to bound the singular 
terms as in (6). 

Let To be the given initial mesh and let S > be a small parameter so that 
#To < S~ d . Later 5 will be chosen such that 5 P « e, where e is the error to be 
achieved between u and uq-, a discrete approximation to u in V7-, and T the 
mesh generated by the algorithm (see proof of theorem 2.1 in section 5.3). Now 
let K £ N be such that 

(K-+l)(27 + d-2) K(2-, + d-2) 

2 2T+d < S < 2 5^+3 . (10) 

Denoting for any element T the elementsize by hr = IT] 1 ^, the constructive 
algorithm reads: 

%),o *~ T) 
3 = 

7, initial (global) refinement to control the error of uq 

7. FIRST LOOP 

do 

Moj = {T£ Tfa : h T > 5} 

<- ref ine(T ^,X j) 
Toj +1 <- complete(7o,j+i) 
3 <- j + 1 
until Mo,j-i = 

J = i 
£ = 1 

7o Selective refinement according to distance to singularities 
7. SECOND LOOP 
while (£ < d(K + 1)) 
&l = \J{T \T £T[ A r T < 2-3} 

2l( T -p-l) , 

= {T C !1{ : h T > 82 <*(2 P+ d) } 
T«+i <- refine (Tf,Mt) 
Tf +1 <— complete(7^ + i) 

end 



The algorithm makes use of two routines that need further explanation. The 
first one, 

Tew <- ref ine(X3i d , M) 

receives a mesh Tid, usually admissible, and a set M. of marked elements from 
Told- It returns a new mesh Tew that is obtained after bisecting once the marked 
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elements according to the newest-vertex bisection rule. The new mesh is not 
necessarily admissible (it may have hanging nodes), but it clearly holds that 

#'?ncw = #^old + #At, 

i.e., #T ncw - #T old = #M. 

The second routine that is used, 

T c *— complete(T) 

receives a mesh T that is not necessarily admissible, and returns a new mesh T c 
which is made admissible by refining the least amount of necessary elements, 
again by the newest-vertex bisection rule. The study of complexity of this 
routine is not as easy as that of the previous one, and it is not true that there 
exists a constant C such that 

< #77 + c(#r £+1 - #77). 

The complexity result that holds — regarding the spreading of refinement implied 
by the completion algorithm — is the following one, which is a little bit weaker, 
but fundamental and sufficient for the purposes of studying optimality of AFEM: 

Theorem 4.1. Let Tq = Tq be an initial admissible mesh of a polygonal (poly- 
hedral) domain Q, in K 2 (M. 3 ), whose elements edges are properly flagged in the 
sense that whenever an interior edge is a refinement edge, it is the common 
refinement edge for all adjacent elements. If the sequence {T^}e>i is obtained 
by subsequent calls to: 

T e+ i <— ref ine(77) 
27+1 <— complete(7^ + i), 

then for k > 1 we have that 

^ i=l ' 
where C is a constant depending only on Tq . 

This result was first proved in [Binev Dahmen DeVore 2004] for triangles, 
and later extended to simplicial meshes of any dimension in [Stevenson 2007] . 

As a consequence of this we have that if now 7oj, T^j Ti, Tf are the meshes 
obtained by our algorithm it holds that 

,d(K+l)-l J-i . 

#t£(k+i) -#%<ci (# T ^+i - # r z c ) + E(# t °j+i - # T °W ■ ( n ) 
^ i=i 3=0 ' 

Remark 4.2. Before proceeding to the proof of our result, some remarks are 
in order: 
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• The idea of the algorithm is to achieve an equidistribution of the error 
following the heuristics stated in the previous section. Since the refine- 
ment is stronger closer to the singularity points, our approach considers 
a sequence of regions Q,£ around them with geometrically decreasing radii 
given by 2~d. The denominator d in the exponent is related to the fact 
that we perform only one bisection to marked elements in refine, and d 
are necessary to achieve a halving of Ht- 

• The algorithm does not take into account the different sizes of the powers 
7i, it just looks at a worst case scenario taking a unified value 7 = m ™' 7l . 
As we will see later, the property 7 > is the only one used in the proof. 
In the same manner, the distance to the singularity points Xi is unified 
by taking the minimum distance symbolized by ry. It may look that 
the simplification introduced by this unification will lead to sub-optimal 
meshes, and it is true that the constant C U ,T in (4) may be bigger than 
necessary with this approach. But this is an a priori approach where we 
want to show the membership of certain functions to the spaces n °t 
caring about the size of their norm. 

• If an efficient construction of the mesh is desired, the algorithm could be 
improved by marking separately according to the different strengths of the 
singularities. This would lead to a better constant C u ,T a , but the overall 
theoretical result will be the same. We decided to present this unified 
approach for the ease of presentation. 

• One important property of the newest- vertex bisection rule is that it leads 
to a sequence of meshes with a uniformly bounded shape-regularity con- 
stant, which depends only on that from the initial mesh 7~ and the new- 
vertex flagging of the initial mesh. We thus have that all the meshes 
T[ obtained by the application of our algorithm are shape-regular with a 
uniform constant. 



4.3 Properties of the Algorithm 

In this section we will bound through a series of lemmas the complexity of the 
resulting mesh T^ K+1 y and in the next section we will relate this complexity to 
the error of the best approximation to u through finite clement functions over 

n-c 

2 d(K+l)- 

The following lemma is related to the termination of the first loop of the 
algorithm in a finite number of steps, and to a control on the number of elements 
added. The termination of the second loop is straightforward, since it is just a 
for loop in disguise. 

Lemma 4.3. The first loop of the algorithm terminates after J iterations, with 
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J < log 2 ^ maXT ^ d r ° + 1 and there exists a constant C\ = 2|f2| suc/i i/ia£; 

(#r 0j+1 - #^.) < d5- d . (12) 



j-i 



7=0 



This implies that for all T £ TJ , |T| < 5 d . 

Proof. Observe first that if one bisects an element T € Tq, J times with J > 

l°g2 ^ maXT ^/° ^ ) + 1, then the measure of the resulting sub-elements will be 

strictly less than S d , and the marking step will not mark them anymore. This 
proves the first part of the statement. 

In order to prove the bound (12) we define, for i > 

T = \t I T £ \Jlo,k A 2M d < \T\ < 2 l+1 (5 d |. 

It is easy to see that even though Ti contains elements belonging to different 
meshes, they do not overlap, and then: 

ioi > J2 \ T \ ^ 6dT = ^ 2i (#^) ; 

which implies that #^ < \^l\b~ d T l . 

Now, applying these estimates, and using that 

oo J J— 1 

(J Ti = {T | T G |J T , fc A \T\ > S d } = (J M 0J) 

i=0 fc=0 j'=0 

we obtain that 

,7-1 J-1 oo 

j=0 7=0 i=0 

and the claim is proved. □ 

Remark 4.4. This proof is a little complicated due to the way the algorithm 
was proposed in order to take into account any previous grading of the mesh. 
Observe that in the first loop we do not refine all the elements, but only those 
which are bigger than the threshold 6, instead of doing just uniform global 
refinements. If we did this, the proof would be simpler, but the number of 
elements in 7~f would be unnecessarily bigger. 

The following lemma is just an observation of the fact that if a point 2 is a 
vertex of a shape-regular triangulation, then the distance of the elements to z 
is an upper bound to the diameter of the element, unless of course the distance 
is zero. This means that the diameter of the elements can grow at most linearly 
with the distance to a point. 
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Lemma 4.5. Let T be a regular mesh such that z is a node, then VT 6 T with 
dist(z,T) ^ we have that \T\ < dist(z,T) d , or h T < dist(z.T). 

This result may be familiar to some practitioners, but it is not completely ob- 
vious. A stronger result was proved in [Nochetto Paolini Verdi 1991, Lemma 5.1], 
but we decided to include its proof here for the sake of completeness. 

Proof. Let T be an element of T and let us define oj t = \J{T | f e T A Tnf ^ 
0}. If z ^ lot then by shape regularity, dist(z,T) > c/it- If z £ ojt\T, then z is 
a vertex of a neighboring element T" and thus dist(z,T) « /ijv ~ /it- D 

The next result implies that the desired grading of the mesh was achieved 
by the algorithm. 

Lemma 4.6. Let T = T^ K+1 y then for < £ < d(K+l) the following property 
holds: 

l , 2€(-,-p-l) 

TeT and r T < 2'* \T\ < S d 2 2 P+ d . 

Proof. We first claim that for each < £ < d(K +1), the following holds for 
the intermediate triangulations Te+i: 

TeT[ +1 and r T < 2 _ * \T\ < <5 d 2 H %Sr Li . (13) 

We prove this by induction on £: By lemma 4.3 it holds for £ = 0. Before 
proceeding, observe that: if T" € T[ and T S 7^, c with k > £: 

TcT' =>■ r T >r T ,. (14) 

Suppose now that (13) holds for ^ and let us prove it for £ + 1. If T G 7^f | _ 2 
and r T < 2 _£ ^ i , there exist T" e 7^ such that T C T', whence r T > < 2~^ , 

2e(~/-p — i) 

and by the inductive assumption \T \ < 8 2 2 v+ d . Now, if already \T'\ < 
Sd 2 ' (e+1 2pl/ 1} t hen the results holds because \T\ < \T'\. Otherwise, T' £ Me+i 
and we have that 

2l(f-p-l) 

, m , l, m ,, 2 2p + d 2«i-fl)( T -p-l) 

|T| < -\T'\ < < S d 2 

because 7 > and d > 2. Thus (13) is proved for £ + 1. 

We now proceed to prove the claim of the lemma: Let TeT such that 
r T < 2-3, then there exist V D T, V £ T[ +l and then by (14), r T > < 2~%, and 

by (13), |T| < \T'\ < ^2 2 "^+d' 1) . □ 

The claim of the previous lemma could have been achieved by simple uniform 
refinement, but this would have destroyed the complexity of the mesh. The next 
lemma shows that the number of marked elements in each iteration is reasonably 
bounded in a way that the overall complexity of the final mesh is under control. 
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Lemma 4.7. There exists a constant C2, depending only on shape regularity, 
such that for 1 < i < d(K + 1); 



#M t = #T £+1 - #T/ < C 2 S- d 2 ip+d _ (15) 

Proof. Recall that in the algorithm wc define Q e = [j{T \ T G T[ A r T < 2~%}, 
and since % + i is obtained from T~[ by refinement only, we have that f2^ = 1J{T 
T € T t+1 : T C whence 



TeT e + 1 , TcQe TeT e+1 \Tf TeT e+1 nT e c TeT l+1 \Tf 

Tcn e 



d 

T- 



But if T € T l+ i\T[, then T is half of an element T" £ A^, and thus by the 
definition of Mi in the algorithm, 



2h d T = h d T , > S d 2 2 P+ d , 

which in turn implies that 

|fi/| > 2 (#T£+1 ~ * T[) = 2 (#Tf+1 ~ #T/) - 

By lemma 4.5 we have that |fig| < C2~ £ and then: 

#T l+1 - #T/ < 2|f^|2^<T d 2~ 2 P+£i < C 2 <5"' i 2 55^R— , 

and the lemma is proved. □ 

The next lemma makes use of the complexity result (11) of the completion 
procedure for the newest-vertex bisection rule, to bound the complexity of the 
final mesh. 

Lemma 4.8. There exists a constant <C 3 , depending only on shape regularity, 
the polynomial degree p, the dimension d, the function u through 7, and Tq, such 
that: 

WI (K+ i)-#To<C 3 S- d . (16) 
Proof. Using (11), lemmas 4.3 and 4.7 we have that 

,d(K+l)-l ,7-1 

^ 1=1 j=0 
,d(K+l)-l 

<CI ^ C 2 S- d 2 — ^ — +Ci5" 



, , -f(2 T + d-2) 

<C5- d C 2 ^2 



-l(2^ + d-2) 

Since 7 > the sum X^=i 2 2 p+ d is finite, and the claim follows taking 

/ -^(2-* + d-2) \ 

C 3 = C C 2 ££Li2 2 P+ d +Cl . □ 
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5 Error 



In this section we bound the best error with finite element functions in terms 
of the complexity of the mesh: 

Theorem 5.1. There exist two constants A.\, A.2, that may depend on u through 
the broken seminorm Wo\ h p+~<- :== \ J2t£T \\D P UaW^tT)) , c%, ki, ji, the 

\\Xi\\wga +1 (n)> WdiWw^in), and the I^g+^-norm of g i; S g <)};,! = l,...,iV, 
the polynomial degree p, the dimension d, shape regularity and T , but otherwise 
independent of K and 6, such that, if T = T^ K+1 y then 

mi Hu-nrll^^A^, (17) 

inf ||u-«r|| lin <A 2 (#T-#T )-§. (18) 

In order to prove this theorem we will consider the regular part uo of u and 
the singular part given by YliLi u i- 

Throughout this section we will use the Lagrange interpolator Iq-Ui of itj, 
which is the finite element function that coincides with Ui at all the nodes, and 
is well defined for each i = 0, 1, . . . , N, since by the assumptions of theorem 4.1, 
all the Ui functions are continuous in SI; see remark 2.2. 

5.1 Estimation of the Regular Part 

Theorem 5.2. There exist two constants C4, C5, depending on the broken 

seminorm \uq\ h p+i^ := (St<eT ^olli^T)) , the polynomial degree p, 

shape regularity and Tq, but otherwise independent of K and 6, such that, if 
T = ^7(K+l)> taen 

\u - Itu \ 1si < C 4 <5 P , 

|«o - Iruo\i, a < c s(# T - #^o)" ? - 

Proof. Since u \t & H P+1 (T) for all T E T , and T was obtained only by 
refinement, u \t G H' P+1 (T) for all T G T, and standard interpolation esti- 
mates [Ciarlet 1978] yield 

K - lTU \l t Q =^2\u - ItUq\\ >t < ^2 h2 T\\ DP+lu a\\ 2 Li{T) 
< S 2p \u \ 2 H ^i {n) , 

where the last inequality is a consequence of lemma 4.3 and the first loop of the 
algorithm. 

Then, by lemma 4.8 

\u - ItuoImi < CiS p < C 5 (#T - #T )-«, 

and the theorem is proved. □ 
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5.2 Estimation of the Singular Part 

Throughout this section, we will denote with u one of the singular terms itj 
defining u in (3), that is, it is defined in polar coordinates around a point x$ in 
f2 as 

u = c l (\n{r l )) k *rfg l (f l )^ (19) 

for some i = 1,2, ... ,N and c,, rj, ki, 7,, <?i, x« as i n t nc assumptions of 
theorem 2.1. 

The three bounds of (6) are the only features of u that will be used in the 
proof of the following theorem. 

Theorem 5.3. There exist two constants Cq, Ct, that depend on the parameters 
defining u in (19), shape regularity and Tq, but otherwise independent of K and 
S, such that, ifT= 1~£, K+1 y then 

\u-I T u\ 1>n <C 6 6 d , 
\u-Itu\ 1jU <C 7 {#T -#T )-$. 

Proof. Let D e = \J{T \ T e T A < dist(x 4 ,T) < 2"^} for < I < 

d(K + 1) and D d[K+l) = \J{T T G T A dist(x*,T) < 2- (K+1 ^}. Then we 
obtain: 

\u - I T u\ 2 i n =^2\u- I T u\ 2 1T 

d(K+l)-l 

= J2 ' u _ J tu|? jT + \u - lTu\ 2 1Dd(K+i) . (20) 

i=0 TcD e 

The second term in (20) can be bounded as follows: 

\u — Itu\ 2 n < Izil? n + \Itu\ 2 n 



x,6T Xi^T 



— : Si + + S3 
From (6) and lemma 4.5, we obtain: 

fc2" 



J 



r2 



-(if+l) 



2^7 / r 27+rf~3 rfr ^ C2 -(K + l)(2 7+ d-2)^ 







For the term B2 we use the fact that on a reference element T, j^ru^ f ^ 
ll IrU IL=°(f) = ll /ru IL°°(T)- B y ( 6 M la; * er,andTc% + i), II Jru IL°°(T) - 
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Ch? T . A proper scaling leads to 



B 2 = £ |/r< T « E h^\(I T u)\ T \l f <C £ 4 7+C 

rcDjiif^i) TC-D<i(jc + i) TC-D<j(jf + i) 

< #{T C D d(Jf+1) : i 4 e T} l^+i)! 22 ^. 



Since for these T's, = 0, lemma 4.6 leads to 

5 2 < #{T C £> rf(x+1) : Xl e T} (2-^(^+i)) 2 ^ < 2 -(^+i)(2 7 +<i-2) ! 

where we have used that the number of elements which have Xj as a vertex is 
bounded by a constant depending only on mesh regularity. 

The term B$ can be bounded using the fact that if dist(xi,T) > 0, then, by 
lemma 4.5, dist(xj,T) ~ |x — Xi\ Vx £ T and thus (6) yields 

|V/ru(x)| < C distix^Ty- 1 < C\x - Xil"/- 1 Vx e T, 

which implies that J T |V/ru| 2 ^ C/ r |a; - Xi| 2 ^ 7-1 ^ dx, and consequently 

B 3 = E / l v/ r«| 2 < C / Ix-Xil^-^dx 



c:2 



< (7 / r 2 (7-i) r d-i ~ c , 2 _(Ar+1)(27+d ~ 2) . 





Combining the three estimates for B\ , B2 and we obtain the following bound 
for the second term of (20): 

|u - I r u\] n < C2- l > K+1 ^+ d -V < CS 2 P+ d . (21) 

Using the usual estimates for the Lagrange interpolator and the fact that u\t G 
H p+1 (T), VT C n\Dd(K-+i) (see remark 2.2), we can bound the first term of (20) 
by: 

d(K+l)-l d{K+l)-l 

E Ei«-h!, t < E E 4 p II^ + HI^ ( t)- ( 22 ) 

fcO TC-Df fcO TC-Df 

Finally, by (6), ifx e T, |£>p +1 u(x)| < C|x-x 4 | 7 ~ p -\ and thus ||DP +1 u|| i2(T)2 < 
Cdist(x t ,T) 2 ^- p - 1 )/i^, by lemma 4.6, h T < 52™^+*? if T e and again 
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by lemma 4.5, we have that 

d(K+l)-l d(K+l)-l 



2 

L 2 (T) 

=0 TCD ( £=0 TcD e 



E Ei-Kt< E E4 p II^ +1 - 



d(X+l)-l 

<C ]T ^ dist^^T) 2 ^" 1 )^ 

£=0 TcD e 
d(K+l)~l 

z c E E 2-^^ 2 Hi2 ^ i2 

£=0 TC-D* 

d(/f+l)-l 

< C<S 2p+d ^ #£>, = CS 2p+d (#T) 

1=0 

Summing up, by (20), (21) and (22), and by lemma 4.8 

I" - *r< n £ ^ 2p+d (#T) = CT 2 P+ d ((#T - #T ) + #T ) 

< cs 2p+d (s- d + #r ) 
<c^<c(#t-#t )-^, 

where we have used that 6 was chosen so that #7o < 6~ d . □ 
5.3 Proof of Main Result 

Proof of Theorem 5.1. Using the estimates of theorems 5.2 and 5.3 we obtain: 
inf ||n-ur|| 1)f2 < inf \u - u r \ i n < C \u - I r u\ i n 



c 



E( u * _ Irui) 



i=0 
N 



1SI 



KC^Km-lTUi^aZCNP, 



i=0 

and then, using lemma 4.8, we have that 



inf \\u-u T \\ in <CN{#T 



□ 



Proof of Theorem 2.1. This is a corollary of theorem 5.1. It is sufficient to 
choose e = Ai5 p . This implies the claim for s small enough, which immediately 
implies the result for all e > 0. □ 
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Remark 5.4. Red-Green refinement. Regarding the other well-known al- 
gorithm for adaptive mesh refinement in two dimensions, namely, the so called 
red-green refinement, the main result presented in this article is still open. How- 
ever, the algorithm stated here can still be used for the construction of the quasi- 
optimal mesh, with obvious modifications due to the fact that a red subdivision 
splits the elements into four sub-elements instead of two. The only remaining 
issue that needs to be solved is to determine if a complexity result bounding the 
spreading of refined elements, similar to theorem 4.1 holds. 
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